Spectra Story

Author

Adam Kemberling

Published

November 5, 2024

Code
# Libraries
library(tidyverse)
library(gmRi)
library(rnaturalearth)
library(scales)
library(sf)
library(gt)
library(patchwork)
library(lmerTest)
library(emmeans)
library(merTools)
library(tidyquant)
library(ggeffects)
library(performance)
library(gtsummary)
library(gt)
library(sizeSpectra)
library(ggdist)
library(pander)
library(ggh4x)
conflicted::conflict_prefer("select", "dplyr")
conflicted::conflict_prefer("filter", "dplyr")
conflicted::conflicts_prefer(lmerTest::lmer)

# Processing functions
source(here::here("Code/R/Processing_Functions.R"))

# ggplot theme
theme_set(
  theme_gmri(
    rect = element_rect(fill = "white", color = NA), 
    strip.text.y = element_text(angle  = 0),
    axis.text.x = element_text(size = 8),
    legend.position = "bottom"))



# vectors for factor levels
area_levels <- c("GoM", "GB", "SNE", "MAB")
area_levels_long <- c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")
area_levels_all <- c("Northeast Shelf", area_levels)
area_levels_long_all <- c("Northeast Shelf", area_levels_long)

# table to join for swapping shorthand for long-hand names
area_df <- data.frame(
  area = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"),
  survey_area = c("SS", "GoM", "GB", "SNE", "MAB", "Northeast Shelf"),
  area_titles = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"))


# Degree symbol
deg_c <- "\u00b0C"



# Function to get the Predictions
# Drop effect fits that are non-significant  ###
get_preds_and_trendsignif <- function(mod_x){
  modx_preds <- as.data.frame(
    # Model predictions
    ggpredict(
      mod_x, 
      terms = ~ est_year * survey_area * season) ) %>% 
    mutate(
      survey_area = factor(group, levels = area_levels),
      season = factor(facet, levels = c("Spring", "Fall")))
  
    # Just survey area and year
    modx_emtrend <- emtrends(
      object = mod_x,
      specs =  ~ survey_area*season,
      var = "est_year") %>% 
      as_tibble() %>% 
      mutate(
        zero = 0,
        non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))
  
    # Preds with signif
    modx_preds %>% left_join(select(modx_emtrend, survey_area, season, non_zero))
  
}
Code
# Broad Distribution
#log2 bins for easy-of-access
#' @title Build Log 2 Bin Structure Dataframe
#' 
#' @description Used to build a dataframe containing equally spaced log2 bins for
#' size spectra analysis. Contains details on the left and right limits, midpoint, bin width, 
#' and a text label for the bins. log2bin number ascends with increasing size for eeasy plotting.
#'
#' @param log2_min 
#' @param log2_max 
#' @param log2_increment 
#'
#' @return
#' @export
#'
#' @examples
define_log2_bins <- function(log2_min = 0, log2_max = 13, log2_increment = 1){
  
  # How many bins
  n_bins  <- length(seq(log2_max, log2_min + log2_increment, by = -log2_increment))
  
  # Build Equally spaced log2 bin df
  log2_bin_structure <- data.frame(
    "log2_bins" = as.character(seq(n_bins, 1, by = -1)),
    "left_lim"  = seq(log2_max - log2_increment, log2_min, by = -log2_increment),
    "right_lim" = seq(log2_max, log2_min + log2_increment, by = -log2_increment)) %>% 
    mutate(
      bin_label    = str_c(round(2^left_lim, 3), " - ", round(2^right_lim, 3), "g"),
      bin_width    = 2^right_lim - 2^left_lim,
      bin_midpoint = (2^right_lim + 2^left_lim) / 2) %>% 
    arrange(left_lim)
  
  return(log2_bin_structure)
}






#' @title Assign Manual log2 Bodymass Bins - By weight
#'
#' @description Manually assign log2 bins based on individual length-weight bodymass 
#' in increments of 1 on the log2 scale. Returns data with bins assigned based on individual
#' length-weight biomass
#' 
#' Uses maximum weight, and its corresponding bin as the limit.
#'
#' @param wmin_grams Catch data prepared for mle calculation, use prep_wmin_wmax
#' @param log2_increment Equally spaced increments to use for log 2 bin sizes. Default = 0.5.
#'
#' @return
#' @export
#'
#' @examples
assign_log2_bins <- function(wmin_grams, log2_increment = 1){
  
  
  #### 1. Set up bodymass bins
  
  # filter missing weights
  size_data <- wmin_grams %>% 
    filter(wmin_g > 0,
           is.na(wmin_g) == FALSE,
           wmax_g > 0,
           is.na(wmax_g) == FALSE)
  
  # Get bodymass on log2() scale
  size_data$log2_weight <- log2(size_data$ind_weight_g)
  
  # Set up the bins - Pull min and max weights from data available
  #min_bin <- floor(min(size_data$log2_weight))
  min_bin <- 0
  max_bin <- ceiling(max(size_data$log2_weight))
  
  
  # Build a bin key, could be used to clean up the incremental assignment or for apply style functions
  log2_bin_structure <- define_log2_bins(
    log2_min = min_bin, 
    log2_max = max_bin, 
    log2_increment = log2_increment)
  
  
  
  # Loop through bins to assign the bin details to the data
  log2_assigned <- log2_bin_structure %>%
    split(.$log2_bins) %>%
    map_dfr(function(log2_bin){
      
      # limits and labels
      l_lim   <- log2_bin$left_lim
      r_lim   <- log2_bin$right_lim
      bin_num <- as.character(log2_bin$log2_bin)
      
      # assign the label to the appropriate bodymasses
      size_data %>% mutate(
        log2_bins = ifelse( between(log2_weight, l_lim, r_lim), bin_num, NA),
        log2_bins = as.character(log2_bins)) %>%
        drop_na(log2_bins)
      
    })
  
  # Join in the size bins
  log2_assigned <- left_join(log2_assigned, log2_bin_structure, by = "log2_bins")
  
  # return the data with the bins assigned
  return(log2_assigned)
  
}

Storycrafting - Northeast Shelf Spectra

I think we need to step back a bit and find the main story again–which is about how size structure of the NES fish community has changed over time, whether that varies by region, and what covariates are associated with those changes.
The Aug draft you developed really focused on comparing regression models to consider drivers based on the fish community size spectrum as viewed from the perspective of the Wigley species set vs. the full species set.
However, I didn’t see figures of those spectra in that draft (maybe there are somewhere else and I’ve lost track). The recent slide deck really focuses on seasonal differences in the size spectra and influence of temperature. - KMills

This document will be split into two sections (the second potentially moving to its own doc).

Part 2: Relationships to External Factors

The aim of this section is to provide context and potentially attribution/correlations that may help explain trends/changes observed in the above section.

I will lead with a characterization of broad regional changes in temperature and fisheries landings totals.

Code
# # Model Dataset for Wigley Species Subset
wigley_bmspectra_model_df <- read_csv(here::here("Data/model_ready/wigley_community_bmspectra_mod.csv"))

# load shelfwide stuff too
wigley_bmspectra_shelf <- read_csv(here::here("Data/processed/shelfwide_wigley_species_bodymass_spectra.csv"))
bot_temps <- read_csv(here::here("Data/processed", "trawl_region_seasonal_bottom_temps.csv"))


# Bolt that on
wigley_bmspectra_model_df <- bind_rows(
  wigley_bmspectra_model_df,
  left_join(wigley_bmspectra_shelf, bot_temps, by  = join_by("est_year" == "year", survey_area, season)))


# Use Wigley bmass spectra model data for plots
wigley_bmspectra_model_df <- wigley_bmspectra_model_df %>% 
  mutate(survey_area = case_when(
    survey_area == "GoM" ~ "Gulf of Maine",
    survey_area == "GB" ~ "Georges Bank",
    survey_area == "SNE" ~ "Southern New England",
    survey_area == "MAB" ~ "Mid-Atlantic Bight",
    survey_area == "Northeast Shelf" ~"Northeast Shelf"),
  survey_area = factor(survey_area, levels = area_levels_long_all),
  season = fct_rev(season))

Bottom Temperature

Seasonal Bottom Temperature Distributions

The Gulf of Maine experiences a much smaller difference in seasonal temperature swings than the other regions. The difference in bottom temperature from Spring to Fall can be smaller than the variability of spring bottom temperature along the whole Northeast Shelf.

The overlap in bottom temperature distributions is very similar to the overlap in mass spectra exponents. Gulf of Maine spectra and bottom temperatures share a good region of overlap, and the regions that experience larger seasonal fluctuations in BT see larger seasonal fluctuations in their size distributions.

Code
# Plot the distribution as a raincloud plot
seasonal_temp_dist <- wigley_bmspectra_model_df %>% 
  ggplot(aes(fct_rev(survey_area), bot_temp, color = season, fill = season)) +
   ggdist::stat_halfeye(
     alpha = 0.4,
     adjust = .5, 
     width = .6, 
    .width = 0, 
    justification = -.2, 
    point_colour = NA) + 
  geom_boxplot(
    fill = "transparent",
    width = .15, 
    outlier.shape = NA) +
  ## add dot plots from {ggdist} package
  ggdist::stat_dots(
    ## orientation to the left
    side = "left", size = 1,
    ## move geom to the left
    justification = 1.12, 
    ## adjust grouping (binning) of observations 
    binwidth = .1) + 
  scale_color_gmri() + 
  scale_fill_gmri() + 
  scale_y_continuous(labels = label_number(suffix = "\u00b0C")) +
  coord_flip() +
  theme(legend.position = "none") +
  labs(y = "Bottom Temperature", x = "", color = "Season",
       title = "Seasonal Bot Temperature")


# Distribution of size distribution exponenents
wigley_b_dist_plot <- wigley_bmspectra_model_df %>% 
  ggplot(aes(fct_rev(survey_area), b, color = season, fill = season)) +
  ggdist::stat_halfeye(
     alpha = 0.4,
     adjust = .5, 
     width = .6, 
    .width = 0, 
    justification = -.2, 
    point_colour = NA) + 
  geom_boxplot(
    fill = "transparent",
    width = .15, 
    outlier.shape = NA) +
  ## add dot plots from {ggdist} package
  ggdist::stat_dots(
    ## orientation to the left
    side = "left", size = 1,
    ## move geom to the left
    justification = 1.12, 
    ## adjust grouping (binning) of observations 
    binwidth = .005) + 
  scale_color_gmri() + 
  scale_fill_gmri() + 
  coord_flip() +
  labs(
    y = "Exponent of Size Distribution", 
    x = "",
    title = "Mass Distributions",
    subtitle = "Bodymass Ranges (min: 1g, max: max observed)")



# Plot next to bottom temperatures
(seasonal_temp_dist |  wigley_b_dist_plot) + plot_layout(guides = "collect")

Seasonal Catchability Differences

The extent that the shape of the size distribution within each region follows a linear relationship to temperature largely is a reflection of the seasonal differences.

If we follow the relationship within a specific season, the relationship is much weaker.

Code
wigley_bmspectra_model_df %>% 
  #filter(survey_area == "Northeast Shelf") %>% 
  ggplot(aes(bot_temp, b)) +
  geom_point(aes(color = season)) +
  geom_smooth(method = "lm", color = "black") +
  facet_wrap(~survey_area, ncol = 1) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(
    title = "Relationship to Bottom Temperature\nReflects Seasonal Catchability Differences",
    subtitle = "Seasonal differences proportional to seasonal temperature range",
    shape = "Survey Area")

Temperature Sentiment

I believe that the size distribution along the shelf is responsive on shorter time scales to temperature changes than hypotheses that predict temperature related changes in growth. This is the seasonal variability in spectra we see.

This suggests to me that the community size distribution is likely less sensistive to long-term trends in temperature since the distribution is able to effectively track swings in temperature that happen within the year.

I think there is an important question about temperature’s role, but I don’t think that local warming is having much impact on the community size distribution.

If it is, it is secondary to other size distribution shifts happening at shorter time scales and across spatial scales.

This raises questions about how seasonal size distribution changes are being achieved.

From the decadal variability work we have reason to believe that individual species are both not tracking temperature changes in space at the rate they are ocurring AND that changes in size at age are happening. These are signals at the species level indicating that temperature effects on growth are ocurring.

There would need to be some compensatory mechanisms present to offset this across the broader community.

Total Landings

The landings information is the newer GARFO data obtained from Andy Beet. It was delivered as containing only finfish, I have removed freshwater species and blue crabs. It should largely reflect marine finfish species at this point.

Code
# From Andy Beet - Crabs removed
landings_annual <- read_csv(here::here("Data/processed/BEET_GARFO_regional_finfish_landings.csv")) %>% 
  select(survey_area,
         est_year = year,
         total_weight_lb = total_live_lb)

shelf_landings_annual <-  mutate(landings_annual, survey_area = "Northeast Shelf") %>% 
    group_by(est_year, survey_area) %>% 
    summarise(total_weight_lb = sum(total_weight_lb, na.rm = T))

# Plot annual totals
bind_rows(list(
  shelf_landings_annual,
  landings_annual)) %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels_long_all)) %>% 
  filter(est_year <= 2019) %>% 
  ggplot(aes(est_year, total_weight_lb/1e6)) +
  geom_col(alpha = 0.5) +
  geom_ma(aes(color = "5-Year Smooth"), size = 1, n = 5, ma_fun = SMA, linetype = 1) +
  scale_color_manual(values = "orange") +
  facet_wrap(~survey_area, ncol = 1, scales = "free") +
  scale_y_continuous(labels = label_number(suffix = "M")) +
  labs(y = "Total Landings (live lb.)", x = "Year", color = "Moving Average")

Zooplankton Community Indices

The ecodata package has changed the zooplankton indices around that it carries. It now contains “regime” metrics for 26 of the main zooplankton taxa.

https://noaa-edab.github.io/tech-doc/zoo_abundance_anom.html

Code
zoo_regimes <- ecodata::zoo_regime
# zoo_anoms <- ecodata::zoo_abundance_anom
zoo_regimes %>% 
  filter(str_detect(Var, "calfin|chaeto|pseudo|oith")) %>% 
  mutate(EPU = factor(EPU, levels = c("GOM", "GB", "MAB"))) %>% 
ggplot(aes(Time, Value)) +
  geom_point(
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(linetype = "5-Year Smooth"),
    size = 1, n = 5, ma_fun = SMA) +
  #scale_color_gmri() +
  facet_grid(EPU~Var) +
  labs()


Driver Correlations

Modeling Driver Relationships

This section is likely to change

The following model has been used to determine the relationship between landings and seasonal bottom temperature on b.

The Northeast Shelf as a whole is not included as its own region in these models, but is plotted below with points just for visual context.

lmer(
  b ~ survey_area*season*scale(roll_temp) + survey_area * scale(log(total_weight_lb)) + (1|est_year), 
  data = wtb_model_df)
Code
#### Model 2: Temperature & Fishing ####


# # Drop NA's
# wtb_model_df <- 
#   drop_na(wigley_bmspectra_model_df, total_weight_lb, bot_temp, b) %>%
#   rename(area = survey_area) %>% 
#   left_join(area_df) %>% 
  
  
wtb_model_df <- drop_na(wigley_bmspectra_model_df, bot_temp, b) %>% 
  # Do rolling BT within a season * region
  group_by(survey_area, season) %>%
  mutate(
    roll_temp = zoo::rollapply(bot_temp, 5, mean, na.rm = T, align = "right",  fill = NA)) %>% 
  ungroup() %>% 
  mutate(
    yr_num = as.numeric(est_year),
    yr_fac = factor(est_year),
    survey_area = factor(survey_area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall")),
    landings = total_weight_lb,
    yr_seas = str_c(season, est_year))



# Make the model
mod2 <- lmer(
  b ~ survey_area*season*scale(roll_temp) + survey_area*scale(log(total_weight_lb)) + (1|est_year), 
  data = wtb_model_df)


# summary(mod2)
# check_model(mod2)

Temperature Relationship

Looking at all the data there is a strong relationship between seasonal size distribution and temperature.

Code
btemp_overall <- ggplot(wtb_model_df, aes(bot_temp, b)) +
  geom_point(aes(shape = survey_area, color = season)) +
  scale_color_gmri() +
  geom_smooth(method = "lm", color = "black") +
  labs(shape = "Area", color = "Season",
       title = "Overall Temperature Relationship")
btemp_overall

But

If we use yearly region * season combinations as our units of comparison the relationship between seasonal bottom temperature and the size distribution falls apart.

Code
# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    mod2, 
    terms = ~ roll_temp*survey_area*season))   %>% 
  mutate(
    survey_area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####
trend_df <- emtrends(
  object = mod2,
  ~survey_area * season,
  var = "roll_temp")


# Drop effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = mod2,
  specs =  ~ survey_area*season,
  var = "roll_temp") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))



# Limit temperature plotting range
actual_range <- wtb_model_df %>% 
  group_by(season, survey_area) %>% 
  summarise(min_temp = min(bot_temp)-2,
            max_temp = max(bot_temp)+2,
            .groups = "drop")



# Plot over observed data
btemp_local <- mod2_preds %>% 
  left_join(actual_range) %>% 
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>% 
  left_join(select(mod2_emtrend, survey_area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = survey_area), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wtb_model_df,
    aes(roll_temp, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(survey_area~., scales = "free") +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Exponent of ISD",
       title = "Within area*season Relationship",
       x = "5-Year Rolling Average Bottom Temperature")
btemp_local

A spring spectra is not any more/less steep if spring is hotter or colder than usual. However, Fall is consistently steeper than spring, and fall is always hotter than spring.

Landings Relationship

Since landings are annual there won’t be any seasonal interaction, only an intercept adjustment for the survey season.

Code
# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    mod2, 
    terms = ~ total_weight_lb * survey_area * season))   %>% 
  mutate(
    survey_area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####
trend_df <- emtrends(
  object = mod2,
  ~survey_area * season,
  var = "total_weight_lb")
#summary(trend_df, infer= c(T,T))


# Drop effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = mod2,
  specs =  ~ survey_area*season,
  var = "log(total_weight_lb)") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T),
    group = str_c(survey_area, "X", season))




# Plot over observed data
mod2_preds %>% 
  left_join(select(mod2_emtrend, survey_area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wtb_model_df,
    aes(total_weight_lb, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(survey_area~., scales = "free") +
  scale_color_gmri() +
  scale_x_continuous(transform = "log10", labels = label_log(base = 10)) +
  labs(y = "Body Mass Spectra Slope (b)",
       title = "Wigley Species, Body Mass Spectra",
       x = "Total Landings (lb.)")

Landings + Temperature Modeling Considerations

With the old landings data we had trouble determining attribution from the models with both temperature and landings simultaneously due to strong collinearities. This may no longer be the case, lets look.

The bigger issue looks to be the relationship between survey season + region and bottom temperature. There are other phenological dynamics (spring plankton bloom, seasonal migrations, spawning, stratification etc.) that happen seasonally that I am hesitant to use temperature alone without controlling for season.

If we check variance inflation factors for the model we currently have nearly all predictors have very high VIF values, beyond what is the typical cutoff.

Which is unfortunate, because it seems like its a pretty reasonable model otherwise.

Code
# Drop effect fits that are non-significant  ###


# summary(mod2)
# summary(temp_landings_mod)
check_model(mod2)

Driver Model Summary

Code
# mod2 %>% 
#   gtsummary::tbl_regression()   %>% 
#   add_global_p(keep = T) %>% 
#   bold_p(t = 0.05)  %>%
#   bold_labels() %>%
#   italicize_levels() %>% 
#   as_gt() %>% 
#   tab_header(title = "Wigley Community - Bodymass Distribution & Drivers")

Side Stories

Everything below is digging into different nuances of the approach or the community composition. These things are of secondary importance and will be moved out.

Spiny Dogfish Migrations

Spiny dogfish constitute a majority fraction of biomass in the MAB during the Spring survey, but are subsequently absent during the fall season.

Their numbers have also been increasing over time and they are unquestionably the dominant “large” species caught in the survey.

Their growing abundance, particularly their spring abundance in the Mid-Atlantic Bight, acts to shallow the spectra as they are one of the larger species routinely sampled by the trawl survey.

If we remove them from the estimation and re-run the slopes this is what changes.

Once removed we see that the Mid-Atlantic Bight is no longer becoming less steep over time, with no long-term trend. Spectra in Georges Bank without the inclusion of spiny dogfish are now declining in both seasons.

Code
# # Run MAB with and without spiny dogfish
# # Run the year, season, area fits for the filtered data
# no_dogs_bmspectra <- group_binspecies_spectra(
#   ss_input = filter(wigley_in, comname != "spiny dogfish"),
#   grouping_vars = c("est_year", "season", "survey_area"),
#   abundance_vals = "numlen_adj",
#   length_vals = "length_cm",
#   use_weight = TRUE,
#   isd_xmin = 1,
#   global_min = TRUE,
#   isd_xmax = NULL,
#   global_max = FALSE,
#   bin_width = 1,
#   vdiff = 2)

# # Save those
# write_csv(
#   no_dogs_bmspectra,
#   here::here("Data/processed/wigley_species_nodogfish_bmspectra.csv"))
Code
# Read them in, do some plotting
no_dogs_bmspectra <- read_csv(
  here::here("Data/processed/wigley_species_nodogfish_bmspectra.csv")) %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels))


# Format a little
no_dogs_bmspectra <- no_dogs_bmspectra %>% 
  mutate(est_year = as.numeric(est_year))


# Join them together
dfish_results <- bind_rows(
   list(
     "Wigley Full" = wigley_bmspectra,
     "Wigley w/o Spiny Dogfish" =  no_dogs_bmspectra), 
   .id = "community") %>% 
  mutate(metric = "bodymass_spectra") %>% 
  mutate(
    survey_area = fct_relevel(survey_area, area_levels))




# Get trends for no dogfish
nodfish_mod <- lmer(
  formula = b ~ est_year * survey_area * season + (1|est_year),
  data = no_dogs_bmspectra)
bmspectra_mod_wig <- lmer(
  formula = b ~ est_year * survey_area * season + (1|est_year),
  data = wigley_bmspectra)

# Put predictions
dogfish_sens_predictions <- bind_rows(list(
  "bodymass_spectra-Wigley Full" = get_preds_and_trendsignif(bmspectra_mod_wig),
  "bodymass_spectra-Wigley w/o Spiny Dogfish" = get_preds_and_trendsignif(nodfish_mod)), 
  .id = "model_id") %>% 
  separate(model_id, into = c("metric", "community"), sep = "-") %>% 
  mutate(metric = fct_rev(metric))



# Plot the differences
ggplot() +
  geom_ribbon(
    data = filter(dogfish_sens_predictions, non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    data = dfish_results,
    aes(est_year, b, color = season, 
        alpha = I(if_else(survey_area %in% c("MAB", "GB"), 1, 0.2))),
    size = 0.8) +
  geom_line(
    data = filter(dogfish_sens_predictions, non_zero == TRUE),
    aes(x, predicted, color = season,
        alpha = I(if_else(survey_area %in% c("MAB", "GB"), 1, 0.2))),
    linewidth = 1) +
  scale_fill_gmri() +
  scale_color_gmri() +
  facet_grid(survey_area ~ metric*community, scales = "free") +
  labs(
    x = "Year",
    y = "Exponent of Size Spectra",
    title = "Shift in Spectra Trend(s) with Exclusion of Spiny Dogfish",
    subtitle = "Spring = Dogfish in MAB, Fall = Dogfish in GB")

Spiny Dogfish Seasonal Dstributions

There is literature documenting a seasonal shift in center of biomass Northward in autumn, and Southward in the spring:

Sagarese et a. 2015.

And another citation documenting their propensity to migrate and range large distances Carlson et al. 20

And some documentation of sexual segregation and habitat usage: Janne B. Haugen et al. 2017

Story Thoughts

Below is a table of some papers operating in the same area or around the same topic that are particularly relevant:

Code
shelf_papers <- tribble(
  ~"Author", ~"Year", ~"Conjecture",
  "Friedland et al.", 2024, "Body size has declined despite more abundance & biomass. Competition is preventing larger sizes.",
  "Krumsick", 2020, "Empirical slopes from trawl survey in labrador sea were steeper than theoretical from nitrogen stable isotope ratios"
)


shelf_papers %>% gt()
Author Year Conjecture
Friedland et al. 2024 Body size has declined despite more abundance & biomass. Competition is preventing larger sizes.
Krumsick 2020 Empirical slopes from trawl survey in labrador sea were steeper than theoretical from nitrogen stable isotope ratios

There are a couple core ideas about what changes have occurred or have been ocurring over the Northeast Shelf with relevance to the community that is sampled by the trawl survey: